subroutine MCrun(nbatch, kinit, lact)
    use util_mc
    use define
    use statistics
    use distlog
    use util_cmfd
    use FMC
    use DF
    use HMGZ
    use util_cmfd_eddington, only : cmfd_EDT
    use nem, only : solve_nem_df, solve_nem_df4
    use mpower
    implicit none
    include "mcbase.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "xsec.FI"
    include "files.FI"
    include "hmgz.FI"
    include "CMFD_LARSEN.FI"
    integer, intent(in) :: nbatch
    real(8), intent(in) :: kinit
    logical, intent(in) :: lact
    integer :: g, i, batch, iter, nsource
    integer, save :: totbat=1
    real(8) :: chk, test, err, shnEtrp, shnEtrp_cmfd
    real(8), save :: wtsum, k_mc, k_cmfd, k_nn, ke, k_trk
    integer, save :: nnt
    integer :: batwp, ntally, ntallywt, nc_super
    logical :: lmwnd
    character*200 logstring
    real(8) :: ts(2), tt(2)
    integer, save :: totiter=0
    
!    integer :: addntrn
    
    ke = 0.0_8
    psifAC = 0.
    nfn = 0.
    
    lmwnd = .false.
    call init_statistics
    
    if (lact .eq. .false.)then
!        if (FDBIACT) bFeedback = .true.
        k_nn = kinit
        k_mc = kinit
        k_cmfd = kinit
!    else
!        if (FDBACT) bFeedback = .true.
    endif
    
    if (swfeedback) then
        logstring = "Feedback on"
        call log_(logstring) 
    else
        logstring = "Feedback off"
        call log_(logstring) 
    endif
    

    !call make_initial_particle(nq, ntrn, bat)
    !call initp_from_dist(nq, ntrn, bat)
#ifdef _2NDHRMN 
    if (lact .eq. .false. .and. swfixtype .eq. 0) then
        psi_f = 1.
        psi_f2(1:totmesh/2) = 1.
        psi_f2(totmesh/2+1:totmesh) = -1.
        
        call add_mpower(psi_f, psi_f2, .true.)
        call comb(history, ntrn(:,bat(NOW)))
        nq(bat(NOW)) = history
    endif
#else
    if (lact .eq. .false. .and. swfixtype .eq. 0) then
        bat(NOW) = NOW
        bat(CUR) = CUR
        bat(NXT) = NXT
        call initp_from_dist(history/swsplitting, ntrn, bat(NOW))
        nq(bat(NOW)) = history
    endif
#endif    
    call cpu_time(ts(1))
    do batch = 1, nbatch
        totiter = totiter + 1
        tempt = 0.
        
	    gk = 0.
	    gktrk = 0.
        absn = 0.
        srcn = 0.
        sctn = 0.
        trk_est= 0.
        jf= 0.
        phi_f2 = 0.
#ifdef _2NDHRMN
        trk_est2 = 0.
#endif
#ifdef CMFDLARSEN
        phi1La=0.
        phi2Ls=0.
#endif
        
        
        ! fmc
        phi_s = 0.
        phi_sa= 0.
        mmt1st= 0.
        mmt2nd= 0.
        mmtAbsr= 0.
        mmtNuF = 0.
                

        if ((NMVWND>0) .and. (totbat .eq.  STMVWND+1)) then
            lmwnd = .true.
            write(*, '(a,i)'), "Movinig window on, size :", NMVWND
            write(hfiter, '(a,i)'), "Movinig window on, size :", NMVWND
        endif 

        if (swfixtype > 0)  then
            selectcase(swfixtype)
                case(1)
                    call initp_from_dist(history/swsplitting, ntrn, bat(NOW))
                    nq(bat(NOW)) = history
                case(2)
                    call initp_from_dist_uniform(history/swsplitting, ntrn, bat(NOW), nq(bat(NOW)))
                case(3)
                    call initp_from_dist_uniform_coarse(history/swsplitting, ntrn, bat(NOW), nq(bat(NOW)))
                case(4)
                    call initp_from_shapeF(history/swsplitting, ntrn, bat(NOW), nq(bat(NOW)))
            endselect
        endif           
        if (delkw>0._8) then
            ke = k_cmfd + shiftk(totiter, delkw)
            k_nn = 1._8/(1._8/k_nn-1._8/ke)
        endif  
        ntally = 0
        nsource = nq(bat(NOW))
        nc_super = 1
        do while(nq(bat(NOW))>0)
            do i=1, nq(bat(NOW))
                call migrate(ntrn(i,bat(NOW)), k_nn*swsplitting, ke*swsplitting, nnt)
            end do
            ntally = ntally + nq(bat(NOW))
            nq(bat(NOW)) = 0
            if (SUPERHST>0)then
                if (nc_super .eq. SUPERHST) then
                    exit
                else
                    call switch_cycle(bat)
                    nc_super  = nc_super +1
                endif
            endif 
            if (ke>0._8 .and. SUPERHST .eq. 0) call switch_wielandt(bat)
#ifdef _WWND
            call switch_wielandt(bat)
#endif 
        end do
        
#ifdef _WWND
        ntallywt = history
#else
        ntallywt = history - nsource + ntally ! weight sum of all neutrons
#endif 
!        call cpu_time(tt(1))

!        do i=1, nq
!            call migrate(ntrn(i,bat), k_nn*swsplitting, ke*swsplitting, nnt)
!        end do
!        ntally = history
!        do while(nqw(batwn)>0)
!            batw = switch(batw)
!            batwn = switch(batwn)
!            ntally = ntally + nqw(batw)
!            do i=1, nqw(batw)
!                call migrate(ntrn_wldt(i,batw), k_nn*swsplitting,  ke*swsplitting, nnt)
!            end do
!            nqw(batw)=0
!        end do

        call leak2netj(ntallywt)
        call evaluate_flux(ntallywt, lact)
        if (lact .eq. .false. .and. RSTDHAT) then
            phird(:,:,batch) = phi_f(:,:)
            jrd(:,:,batch) = jf(:,:)
        endif
#ifndef _PUREMC
#ifdef _2NDHRMN 
        call add_mpower(psi_f, psi_f2, .false.)
        call regen_flux(phi_f, phi_f2nd)
#endif
#endif

        !if (FMSHAPE) CALL update_fmshape(psifAC)
        !if (FMSHAPE) CALL update_fmshape(psi_f)
        
        !if (lact) call write_distlog1
#ifndef _NOLOG
        call write_distlog1
#endif
!        write(*,*) sum(phi_f(:,1))
        
!        jf = jf*wtsum
!       ! check nodal balance : fixed source balance
!       
!        do g=1, ng
!            do i=1, totmesh
!                chk = jf(2*i-1, g) +jf(2*i, g) + absn(i, g) - srcn(i, g) + sum(sctn(g, :, i))-sum(sctn(:, g, i))
!                if (abs(chk)>1e-10) then
!                    test = 0.
!                end if 
!            end do
!        end do
!        jf = jf/wtsum
      
        call estimateK(ntallywt, k_mc, k_trk)
        shnEtrp = ShannonEntropy()
        if (ntally .ne. nsource) then
            !k_mc = (k_mc*ke)/(k_mc+ke)
            write(*, '(i6,2x,i10,2x,i10,2x,f7.5,2x,f7.5,2x)',advance="no"),batch, nsource, ntally, k_mc, k_trk
            write(hfiter, '(i6,2x,i10,2x,i10,2x,f7.5,2x,f7.5,2x)',advance="no"),batch,  nsource, ntally, k_mc, k_trk
        else
            write(*, '(i6,2x,i10,2x,f7.5,2x,f7.5,2x)',advance="no"),batch, nsource, k_mc, k_trk
            write(hfiter, '(i6,2x,i10,2x,f7.5,2x,f7.5,2x)',advance="no"),batch,  nsource, k_mc, k_trk
        end if
#ifndef _PUREMC
#ifndef _2NDHRMN 
        k_cmfd = k_mc
        if (NMVWND .eq. 0 ) then
!            if (lact .and. RSTDHAT) then
!                call homogenize_acc_reduce(batch, inactive)
!            else
                call homogenize_acc
!            endif
            call sethk(avgk_mc)
        elseif (NMVWND .eq. 1) then
            call homogenize_cycle
            call sethk(k_mc)
        else
            call homogenize_wnd
        endif 
        
        call st_avgk(k_mc)
!        call balance_control
        ! cmfd step
!        call cmfd_EDT(k_cmfd, err, iter, lmwnd)
        call cmfd(k_cmfd, err, iter, lmwnd, lact)
        if (NMVWND .eq. 1) then
            fphi1_hf=0;gphi0_hf=0;phi0=0;phi2=0;phi2_hf=0;phi1p=0;phi1c=0;fphi2fs=0;phi2f=0;
        endif
!        call add_statistics(k_mc, k_cmfd)
        call st_avgkcmfd(k_cmfd)
        
        
!        call cpu_time(tt(2))
!        gtally(1) = tt(2)-tt(1)
!        tt(1)= tt(2)
#ifndef _NOFMC
        ! FMC
        call slove_fmc_tdiag(lact)
        
!        call solve_nem_df
!        call solve_nem_df4
#endif        
#endif
#endif

!        call evaluate_corcoef(batch, k_mc, phi_f(mpeak,1), gtally_temp)
        ! if (lact) call write_distlog2
        call write_distlog2
        shnEtrp_cmfd = SHN_CMFD()
!        gtally(3)=err_ref
!        write(*,*) sum(phi_f(:,1))
#ifdef _NOACC
        write(*, '(6(f12.5,2x))'), avgk_mc, stdk_mcc, stdk_mc, avgk_cmfd, stdk_cmfd, shnEtrp
        write(hfiter, '(6(f12.5,2x))'), avgk_mc, stdk_mcc, stdk_mc, avgk_cmfd, stdk_cmfd, shnEtrp
#else
        call cpu_time(ts(2))
#ifndef _TRPCHK
        write(*, '(2(f7.5,2x),f7.2, 2x, f7.5, 2x, f7.5, es12.3, es13.5, es13.5, es13.5)'), avgk_mc, stdk_mcc, stdk_mc*1e5, shnEtrp, shnEtrp_cmfd, ts(2)-ts(1), gtally(1), gtally(2), gtally(3)
        write(hfiter, '(2(f7.5,2x),f7.2, 2x, f7.5, 2x, f7.5, es12.3, es13.5, es13.5, es13.5)'), avgk_mc, stdk_mcc, stdk_mc*1e5, shnEtrp, shnEtrp_cmfd, ts(2)-ts(1), gtally(1), gtally(2), gtally(3)
#else
        write(*, '(2(f7.5,2x),f7.2, 2x, f7.5, e12.3, e12.3, e12.3, e12.3, e12.3, e12.3, e12.3)'), avgk_mc, stdk_mcc, stdk_mc*1e5, shnEtrp, ts(2)-ts(1), gtally(1), gtally(2), gtally(3), tempt(1), tempt(2), tempt(3)
        write(hfiter, '(2(f7.5,2x),f7.2, 2x, f7.5, e12.3, e12.3, e12.3, e12.3, e12.3, e12.3, e12.3)'), avgk_mc, stdk_mcc, stdk_mc*1e5, shnEtrp, ts(2)-ts(1), gtally(1), gtally(2), gtally(3), tempt(1), tempt(2), tempt(3)
#endif        
        ts(1) = ts(2)
        
#endif
 
       
!        call cpu_time(tt(2))
!        gtally(2) = tt(2)-tt(1)
!        tt(1)= tt(2)
#ifdef _2NDHRMN 
        call comb(history, ntrn(:, bat(NXT)))
        nq(bat(NXT)) = history
#else
        !** process for next step **!
        if (swfeedback) then
        if(lact .eq. .true. .or. batch > delay_cmfd) then
#ifdef _FMCFDB
#ifndef _NOFMC
            call feedback_FMC(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
#else
#ifdef _RSMP
            call feedbackRSMP(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
#else
            call feedback(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
#endif
#endif

#else
#ifdef _RSMP
            call feedbackRSMP(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
#else
            call feedback(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
#endif
#endif
            if (swfixtype .eq. 4) call update_pdf_coarse(psi_c)
        else
            call init_weight(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
        end if
        else
            call init_weight(nq(bat(NXT)), history, ntrn(:, bat(NXT)))
        end if
                
!        call cpu_time(tt(2))
!        gtally(3) = tt(2)-tt(1)
!        tt(1)= tt(2)

        totbat = totbat+1
#endif

!        k_nn = k_mc
        if ((swfeedback .eq. .false.) .and. (swfixtype .eq. 0)) then
!            k_nn = k_mc*nqn/history
            k_nn = k_mc
!             k_nn = k_mc
        else 
!            k_nn = k_mc
            k_nn = k_mc
        endif
        
!        call solve_DF
        
        call switch_cycle(bat)
        
        if (MOREPTC) then
            if (lact .eq. .false. .and. batch .eq. nbatch-1) then
                call adjustWeight(wtsum, history, ntrn(:,bat(NXT)), nq(bat(NXT)), NMPTCR)
            endif
            if (lact .and. batch .eq. NMCYCLE-1) then
                call adjustWeight(wtsum, history, ntrn(:,bat(NXT)), nq(bat(NXT)), 1._8/NMPTCR)
            endif        
        END IF 
        if (WRTREF) then
            open(2001, file='REF.txt', form='binary', status='unknown')
            write(2001) phifAC
            close(2001)
        endif 
    end do
    
end subroutine


subroutine MC(initk)
    use allocs
    use define
    use distlog
    use util_mc
    use mvwnd
    use FMC
    use DF
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    include "files.FI"
   
    integer :: nbatch, seed, ichar1st, i
    real(8), intent(in) ::initk
    real(8) :: tchk(2), ct, dummy
    character*200 logstring, strseed
!    integer :: i, g, i_c

    ! FMC
    call init_fmc
    
    ! DF
    call init_df
    
#ifdef RNGMCNP
    call initRN()
    logstring = "RNG : MCNP"
    call log_(logstring) 
#endif    
#ifdef RNGMCCARD
    logstring = "RNG : MCCARD"
    call log_(logstring) 
#endif
    if (impl_cpt) then
        logstring = "Implicit capture on"
        call log_(logstring) 
    endif
#ifdef _EXPTR
        logstring = "Exponential transform on"
        call log_(logstring) 
#endif    

    call getarg(1, strseed)
    ichar1st=ichar(strseed(1:1))
    read(strseed, '(i)'), seed
    call initRN(seed)
!    do i=1, seed
!        dummy = GetRN()
!    end do

    ! initialize moving window 
    if (NMVWND>0) call init_mvwnd(NMVWND, totmesh, totm_cmfd*2, ng)
    
#ifdef _LANISO
    logstring = "Anisotropic scattering with Larsen's Xsection"
    call log_(logstring)
    call init_anisoSct
#endif


    
    call init_distlog
    call init_util_mc
    if (FMSHAPE) call update_fmshape(psi_f)
    
    if( PDFREF ) phi_f = refphi
    if( FLTIGS ) phi_f = 1.

    call init_pdf
    call make_pdf(phi_f)
    if ( FXDSRC .ge. 3) then
        call init_pdf_coarse
        call make_pdf_coarse(phi_f)
    endif
    if (FXDSRC .eq. 4) call formfunction(phi_f)
        

    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! inactive cycle 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call prp_inactive
    call cpu_time(tchk(1))
    call MCrun(inactive, initk, .false.)
    call cpu_time(tchk(2))
    ct = tchk(2)-tchk(1)
    write(logstring, *) "time(sec) : ", ct
    call log_(logstring)
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! active cycle 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call prp_active
    call cpu_time(tchk(1))
    call MCrun(active, initk, .true.)
    call cpu_time(tchk(2))
    ct = tchk(2)-tchk(1)
    write(logstring, *) "time(sec) : ", ct
    call log_(logstring)
    call finalize_distlog
    call write_result
    
    call check_bias
    
end subroutine


subroutine prp_inactive
    use define
    implicit none
    include "pspec.FI"
    include "mcbase.FI"
    include "arrays.FI"
    include "accm.FI"
    character*200 logstring
    
    bat(NOW) = NOW
    bat(CUR) = CUR
    bat(NXT) = NXT

    jfAC     = 0.
    pjAC     = 0.
    psifAC   = 0.
    phifAC   = 0.
    phif2AC  = 0.
    phisAC   = 0.   
    phisaAC  = 0.  
    mmt1stAC = 0. 
    mmt2ndAC = 0. 
    mmtAbsrAC = 0. 
    mmtNuFAC  = 0. 
!    trk_estT = 0.

    swfeedback = .false. 
    if (FDBIACT) swfeedback = .true.
    
#ifdef CMFDTEST
    swfixtype = 0
    swfmshape = 0
    swsplitting = 1
#else
    swfixtype = FXDSRC
    swfmshape = FMSHAPE
    swsplitting = SRCSPLT
#endif

    logstring = "Inactive cycle"
    
    call log_(logstring)
    if (swfixtype > 0) then
        write(logstring, '(a,i01)'), "Fixed Source MC in Inactive Cycle  type:",swfixtype 
        call log_(logstring)    
    endif
    
    if (swsplitting .eq. 2) then
        logstring = "Source Spliting"
        call log_(logstring)    
    endif
    if (swfmshape) then
        logstring = "Linear Approximation of flux in a fine mesh"
        call log_(logstring)    
    endif
end subroutine


subroutine prp_active
    use util_MC, only : write_result
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "stsol.FI"
    include "CMFD_LARSEN.FI"
    character*200 logstring
    swfeedback = .false. 
    if (FDBACT) swfeedback = .true.
#ifndef CMFDTEST
    swfixtype = 0
    swfmshape = 0
    swsplitting = 1
#else
    swfixtype = FXDSRC
    swfmshape = FMSHAPE
    swsplitting = SRCSPLT
#endif
    
    if (RSTDHAT) then
        logstring = "Reset Dhat."
        call log_(logstring)
        phi1Lah=0
        phi2Lsh=0
        phi1Laacc=0
        htjr=0
        jfAC     = 0.
        pjAC     = 0.
        psifAC   = 0.
        phifAC   = 0.
        phif2AC  = 0.
        phisAC   = 0.   
        phisaAC  = 0.  
        mmt1stAC = 0. 
        mmt2ndAC = 0. 
        mmtAbsrAC = 0. 
        mmtNuFAC  = 0. 
        fphi1_hf=0;gphi0_hf=0;phi0=0;phi2=0;phi2_hf=0;phi1p=0;phi1c=0;fphi2fs=0;phi2f=0;
    endif
    logstring = "Active cycle"
    
    call log_(logstring)
    if (swfixtype > 0) then
        write(logstring, '(a,i01)'), "Fixed Source MC in Inactive Cycle  type:",swfixtype 
        call log_(logstring)    
    endif
    
    if (swsplitting .eq. 2) then
        logstring = "Source Spliting"
        call log_(logstring)    
    endif
    if (swfmshape) then
        logstring = "Linear Approximation of flux in a fine mesh"
        call log_(logstring)    
    endif
    stmc=0
    stsqmc=0
    
    xsqsum=0;xysum=0;xsum=0;ysum=0
end subroutine

subroutine check_bias
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer :: i, g
    real(8) :: xbar, ybar, xsqexp, xyexp
    real(8) :: ract
    ract=1._8/active
    do g=1, ng
        do i=1, totm_cmfd
            xbar=xsum(i,g)*ract
            xsqexp=xsqsum(i,g)*ract
            ! for the left
            ybar=ysum(i*2-1,g)*ract
            xyexp=xysum(i*2-1,g)*ract
            beta(i*2-1,g)=xsqexp/(xbar**2)-xyexp/(xbar*ybar)
            
            ! for the right
            ybar=ysum(i*2,g)*ract
            xyexp=xysum(i*2,g)*ract
            beta(i*2,g)=xsqexp/(xbar**2)-xyexp/(xbar*ybar)
        enddo
    enddo
    
    open(990,file ="bias.txt",status='unknown')
    do g=1, ng
        write(990, '(2000(e15.8,1x))') beta(1:totm_cmfd*2,g)
    enddo
    close(990)
    
endsubroutine